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Abstract This paper proposes a probabilistic approach for the detection and 
the tracking of particles in fluorescent time-lapse imaging. In presence of very 
noised and poor quality data, particles and trajectories can be characterized by 
ana-contrario model, that estimates the probability of observing the structures of 
interest in random data. This approach, first introduced in the modeling of human 
visual perception and then successfully applied in many image processing tasks, 
leads to algorithms that do not require a previous learning stage, nor a tedious 
parameter tuning and are very robust to noise. Comparative evaluations against a 
well established baseline show that the proposed approach outperforms the state 
of the art. 

Keywords particle detection • particle tracking • a-contrario approach • 
time-lapse fluorescence imaging 


1 Introduction 

Advances in microscopy and fluorescence technology over the last years, have led 
to the collection of huge amounts of fluorescent biological data, which require 
the use of automatic image processing tools to be analyzed quantitatively [lEl 
El El ESI- However, mainly due to the poor data quality and the complexity 
of subcellular component dynamics, fluorescent time-lapse imaging represents a 
challenging domain for automatic image analysis. Indeed, when a single point- 
source of light is brought to a focus with a lens, the point-source image has a typical 
normalized intensity distribution called Point Spread Function (PSF), which can 
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Fig. 1: Example of particles corresponding to protein aggregates into 
bacteria cells visualized through a fluorescence microscopy. 


be very accurately approximated by a 2D Gaussian m- In addition, in fluorescent 
imaging, the point-sources are often subcellular structures, whose size is typically 
smaller than the resolution limit of the microscope, resulting in a diffraction limited 
spatial resolution. As a consequence, the objects of interest typically appear as 
bright spots severely blurred, commonly called particles (see Figj^ over a possible 
widely varying background intensities. Furthermore, specially in live cell imaging, 
the signal-to-noise ratio (SNR) is typically very low. Indeed, the intrinsic photon 
noise introduced in the imaging process can be reduced only by increasing the light 
intensity or the exposure time, which in turn causes the fading of the fluorescent 
signal, a process called photobleaching m- To the poor data quality, the analysis 
of time-lapse sequences adds further challenges since it often requires to track 
over time multiple appearing/disappearing particles which undergo heterogeneous 
motion. In addition, tracking algorithms have to cope with missed detections and 
spurious particles that may arise from the detection step. 

At present, a large number of algorithms and software tools are available for the 
spatial detection and temporal linking of particle and, in recent years, a number 
of works have attempted to address the problem of objectively assessing available 
methods under different experimental conditions [ISIIIIIEIIISIIIIT], but they were 
limited to either one aspect of the task (particle detection or particle temporal link¬ 
ing) or to a single biological scenario. A recent work published on Nature Methods 
by Chenouard et al. m has collected the results of an open competition organized 
in 2012 to which participated 14 teams. The challenge was to test tracking algo¬ 
rithms of each team on common simulated data sets, representative of different 
biological scenarios, and to evaluate their performances by using a common set of 
evaluation criteria. This study indicates that, at present, there exists no univer¬ 
sally best method for particle detection and tracking since a method reported to 
work for certain experiments may not be the right choice for another application. 
In particular, it has been shown that most available tracking techniques cannot 
cope with high levels of noise and high particle density. 

In this manuscript, we propose a probabilistic approach for the detection and 
the temporal linking of near-circular particles in biological images and we show 
its advantages over existing ones in terms of control of false alarms, robustness to 
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noise, and reduced number of tuning parameters. In addition, being the approach 
unsupervised, it avoids the drawback of supervised methods such as the bias- 
variance dilemma and the need of a cumbersome learning stage. The contribution 
of this paper is twofold: first, we propose a novel method for particle detection 
based on the a-contrario framework ED; second, using the baseline issue of the 2012 
Particle Tracking Challenge, we evaluate the performances of a particle tracking 
method that takes the proposed method for particle detection as input and links 
particle in successive frames by a recently introduced method for particle temporal 
linking [T], also based on the a-contrario framework. 

In the next section, we introduce the state of the art on particle detection and 
particle temporal linking. In sectionwe recall the formalization of the a-contrario 
framework that will be used in section |^and section [^to explicit the a-contrario 
model for particle detection and particle tracking respectively. We devoted section 
I^to the introduction and discussion of the experimental results. Finally, in section 
^ we draw our conclusions. 


2 Related work 

2.1 Particle detection 

Particle detection methods can be broadly classified into supervised and unsu¬ 
pervised. Basically, supervised methods learn the particle model appearance from 
annotated training data consisting of positive and negative samples, whereas un¬ 
supervised methods assume some particle appearance model and rely on different 
filtering and detection techniques. Typically, in unsupervised methods, the parti¬ 
cle model derives from Gaussian approximations of the PSF m niiMi ESI E2] 
or from wavelet decompositions musKis], or from feature-based approaches m 
or from mathematical morphology [SHlIsaiMlISQ]. 

Methods deriving from a Gaussian approximation of the PSF include the Top- 
Hat Filter (TH) [TOl [9] and the Spot-Enhancing Filter (SEF) [53]. TH [TOl 0 
extracts small, compact or rounded objects from images. This is achieved by ex¬ 
ploiting apriori information about object shape and predetermined information 
about their intensity from a circular interior region around a candidate point and 
a surrounding annular region. If the brightness difference in the two regions ex¬ 
ceeds a threshold level, the candidate point is considered to be a particle. SEE 
was proposed by Sage et al. [53] and basically consists of an enhancement filter 
resting on estimation theory, following which the maximum SNR detector of a 
given signal, or template, in additive stochastic noise is provided by the whitened 
matched filter. The matched filter is obtained by convolving the unknown signal 
with a conjugated time-reversed version of the template. The authors showed that 
the optimal detector of Gaussian-like particles in a fractal-like noise with a spec¬ 
tral power density that decays like l/cj^, where uj is the radial spatial frequency, 
is the Laplacian of a Gaussian, also known as Mexican hat filter. 

Methods deriving from wavelet decomposition include the Wavelet Multiscale 
Product (WMP) and the Multiscale Variance-Stabilizing Transform Detector (MSVTV). 
WMP is based on the wavelet decomposition introduced in m. under which real 
objects, contrary to noise and randomly distributed data, are represented by a 
small number of wavelet coefficients that are correlated and propagated across 



4 


Mariella Dimiccoli et al. 


scales. Hence, objects can be detected simply by thresholding the multiscale prod¬ 
uct image. MSVTV is based on the multiscale variance-stabilizing transformation 
m- Since objects under this transformation have to be localized in both space 
and frequency, large values of the transformed image usually correspond to some 
structure and smaller ones to noise. The significant wavelet coefficients are de¬ 
tected by performing multiple hypothesis testing through the Benjamin-Hochberg 
procedure izaj. In the reconstructed image, only significant coefficients are nonzero 
and the background is largely removed whereas objects are preserved. Detection 
is performed by thresholding the reconstructed image. 


Methods deriving from morphological operators include the Grayscale Open¬ 
ing Top-Hat Filter (MTH) and H-Dome Based Detection (HD). The MTH [5811^ 
is obtained by subtracting to the image its opened version, obtained by using a 
fiat disk as structuring element. The subtraction yields an image with only the 
removed objects which correspond to round light objects on a dark background. 
Contrary to the MTH, which select only compact structures smaller that the struc¬ 
tural element, HD [SSIEO] acts by subtracting from the original image / the mor¬ 
phological reconstruction of the image / — h, where h is a constant image. This 
detector depends on the local contrast regardless the morphology or the scale of 
the objects. A threshold above h then only keeps the h-contrasted peaks (i.e lo¬ 
cal maxima). Its shortcoming is that small contiguous particles are extracted as 
one connected region because the size of the structuring element is wider than 
the minimum distance between the peaks of adjacent particles. Methods deriving 
from feature-based approaches are introduced and called Image Features Based 
Detection (IDF) in [61]. The key idea underlying these methods is to combine 
image intensities with local curvature information. This is achieved by computing 
at each pixel the determinant of the Hessian matrix with a smoothing scale m 
or, alternatively, by multiplying the value of the determinant of the Hessian with 
the intensity values. Supervised methods for particle detection include AdaBoost 
(AB) and Feature Discriminant Analysis (FDA). Both approaches work by 
classifying image patches extracted from the images through a sliding window ap¬ 
proach as particle or background. A set of four Haar-like features are extracted 
from each image patch. The AB classifier consists of a sequence of weak classifiers 
which are combined in a weighted sum to create the final output of the boosted 
classifier. FDA [33] is a statistical technique that aims at finding, during training, 
a projection where the class separation (particle and background) is maximized, 
taking into account the mean and the covariance matrix for each class. This in¬ 
formation is used during testing to generate a classification map, which convey 
particles when thresholded at a value, which is also automatically estimated from 
training data. From the outstanding quantitative comparison work made by Smal 
et al. [61], which includes seven unsupervised (TH, SEF, WMF, MSVTV, MTH, 
HD and IDF) and two supervised methods (AB and FDA), differences in perfor¬ 
mances are negligible at high SNRs (> 5) but performances of most methods drop 
out at SNR lower than 4. Taking into account also the number of parameters and 
the sensitivity of the methods to parameter changes, supervised methods achieve 
overall better performances at low SRNs but at the price of a cumbersome training 
stage, which may possibly introduce a bias. 
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2.2 Particle temporal linking 

Early methods for particle temporal linking have addressed the easier problem of 
tracking a single-particle but they may also track multiple particles whose trajec¬ 
tories are sufficiently separated in space [HminiisTiiiEZ]. Existing algorithms for 
tracking multiples particles can be broadly classified into deterministic and prob¬ 
abilistic approaches. Typically, deterministic approaches m usi ini ill ES] act 
by first estimating the position of each particle in each frame independently and 
then by linking particles in successive frames. Instead, probabilistic approaches 
[n SOI ESI mum ED include a spatio-temporal filtering mechanism which allows 
to better exploit temporal information capturing the uncertainty of the measures 
due to noise (random variations) and other inaccuracies. Deterministic approaches 
may be local or global. Local approaches to link particles are willing to fail when 
particles move quickly, close to each other and in presence of spurious/missed 
detections. Deterministic global strategies such as the global nearest neighbor ap¬ 
proach [221 sam] or Gaussian template matching m attempt to find and to 
propagate the single most likely hypothesis at each frame. In multiple hypothesis 
tracking m, temporal information is exploited to solve assignment ambiguities by 
delaying the association task between a set of measurements and a set of tracks to 
future observations that will resolve the conflict. This approach suffers from com¬ 
binatorial explosion. Typically, deterministic global strategies assume a motion 
mode of the particles and therefore are enable to deal with a variety of trajectory 
speeds at the same time. A deterministic global approach that does not assume 
motion modes was proposed by Sbalzarini and Koumoutsakos [55]: associations 
between points corresponding to the same physical particle in subsequent frames 
are computed by minimizing a cost functional with topological constraints through 
a greedy algorithm. Another example of global deterministic approach was pro¬ 
posed by Jaqaman et al. m- This method works by first linking particles between 
consecutive frames and then by linking the resulting track segments into complete 
trajectories. Both steps are formulated as global combinatorial optimization prob¬ 
lems whose solution identifies the overall most likely set of particle trajectories 
throughout the sequence. 

More sophisticated probabilistic multi-particle tracking algorithms model the 
object trajectory as a dynamic system, whose state-space evolution over time is 
described through two equations: the measurement equation relating the observed 
data Z to the state vector X and the system transition equation for the state 
vector X. The state vector X represents the object motion to be estimated and 
the measurement vector Z represents the observed motion [SIIMIIIS]. A popular 
probabilistic approach is the Bayesian Sequential Estimation, which allows the re¬ 
cursive estimation of the so called filtering distribution p{Xt\Zi:t), describing the 
object state conditional on all the observations seen so far. If the posterior density 
at every time step is Gaussian and the system model along with the measurement 
model are linear, then the estimation can be done optimally by a Kalman filter [5]. 
In contrast to Kalman Eilter, Particle Eilter [SSIIIS] exploits better the temporal 
information and allows to approximate models that are not linear and/or not Gaus¬ 
sian. The filtering distribution is presented as a set of samples, or particles, with 
associated weights. The weights are propagated over time to give approximations 
of the filtering distribution at subsequent time steps. An important shortcoming 
of Particle Eilters is that it leads to spurious detections. To address this problem. 
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Godinez et al. m proposed to a probabilistic data association approach that com¬ 
bines a top-down strategy driven by the Kalman filter and a bottom-up strategy 
using standard localization algorithms for fluorescent particles. Other probabilistic 
approaches include probabilistic variations of multiple hypothesis tracking (MHT) 
[671 EH [Hilo], multitemporal association tracking m and probabilistic data as¬ 
sociation m- In probabilistic MHT m, instead of assigning measurements to 
tracks as in traditional MHT algorithms, the probability that each measurement 
belongs to each track is estimated using a maximum a posteriori method. This 
algorithm has poor performances in cluttered environments. To solve this problem 
Chenouard et al. m proposed a MHT approach where detections are linked to 
form target trajectories by using a Bayesian framework aiming at building the set 
of tracks that maximizes the likelihood of the associations between tracks and mea¬ 
surements from the images of the sequence. In multitemporal association tracking 
[6^ the tracking problem is posed in a graph-theoretic framework where the de¬ 
tections are treated as vertices of a graph and the edges are possible inter-frame 
associations. Instead of finding associations between detected objects, the associa¬ 
tion is found between feasible paths and the current set of tracks. The association 
is done by minimizing a cost function that approximates the Bayesian a posteri¬ 
ori probability estimate for the data association problem. The key difference with 
other approaches to solving the multitarget tracking problem is that it allows a 
single track to be assigned to more than one path. Shafique et al. m perform 
probabilistc data association between set of particles belonging to different frames 
by looking for the maximum matching of a bipartite graph, where the partite 
set correspond to the set of points detected in two successive frames. The greedy 
algorithm has the advantage of allowing the use of different motion models and 
cost functions. A major drawbacks of these methods are the assumptions about 
the probability distributions that do not necessarily hold and the large number of 
parameters. 

A comparative evaluation of virus tracking methods has been done by Godinez 
et al. m- They provided a performance evaluation of eight approaches suggest¬ 
ing that probabilistic approaches yield better performances than deterministic 
approaches. 


3 The a-contrario framework 

The a-contrario framework rests on a perception principle stated by Helmholtz 
following which the human visual system detects structures in a group of objects 
when their configuration, according to one or several Gestalt laws, are very unlikely 
to happen by chance in a random setting. Basically, the a-contrario methodology 
requires two ingredients: a naive model, that describes typical situations where 
no structure should be detected and one or several measurements defined on the 
structures of interest. For instance, when trying to discover alignments of points 
in an image, the naive model should consist in a uniformly and independent draw 
of points where no alignments should be detected (see Fig. |^. The measurements 
should define in what way an observation can be significant and are usually related 
to the visual saliency of the structure. For instance, for the image in Fig. the 
alignments of points should pop out because, assuming the naive model, they are 
not likely to happen by chance considering the total number of points in the image. 
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More formally, if the measurement function is high when the structure is pregnant, 
the amount of surprise when observing the measurement x can be related to the 
probability P(X > x), where X is the random variable corresponding to the 
distribution of x in the naive model. We will usually have several measurements 
and in the classicala-contrario framework the amount of surprise will be measured 
by a Number of False Alarms (NFA), defined formally as follows. 

Definition 1 (Number of False Alarms) Let (W)i<z<Ar he a set of random 
variables. A family of funetions [Fi h))i is a NFA (number of false alarms) for 
the random variables {Xi)i if 

V6 > 0, E(#{f, Fi(V) < e}) < 6 (1) 

(as usual, the notation stands for the eardinal of the set S). 

The NFA ensures that the average number of detections made in the naive model 
(false detections) at level e is less than e. Such detections are said meaningful. 



Fig. 2: Illustration of Hemholtz principle. According to Helmholtz principle, 
we a priori assume that the dots should have been drawn from an uniform distri¬ 
bution as in the right image. However, in left image, we perceive a group of aligned 
dots because such a structure is very unlikely to happen by chance in a random 
setting. Actually, also in the right image there is an alignment of three points but 
it does not pop out, because it is likely to happen by chance considering the total 
number of points. 


4 A-contrario particle detection 

To apply the a-contrario framework to the detection of particles, we need to spec¬ 
ify the naive model Ho as well as a statistical measurement function m able to 
characterize the visual saliency of the particles we are looking for. As naive model, 
we take the realization of a Gaussian stochastic process with mean /i and standard 
deviation a {Ho = A/’(yn; a)). In such random image, all image pixels are indepen¬ 
dent random variables with uniform probability distribution over some interval 
and therefore no structure should be detected. As measurement function to be 
performed at any given location (x, y) of the image grid T , we consider the local 
contrast of a small patch centered at (x, y) with respect to its local background. 
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As observed by Grossjean and Moisan m, modeling the local context is neces¬ 
sary to bypass the sensitivity of the observer to low frequencies in the detection 
task. To this goal, we define the statistical measure m as the difference of one 
principal measurement, say mi, defined on a disc of radious R centered at (x,y) 
representing the detection area, and a context measurement, say m 2 , defined on a 
ring surrounding the previous disc of radious R, with (a > 1) (see Fig. [^. Since 
a particle should be characterized by an high measurement function in the inner 
disc and a low measurement function in the outer disc, assuming that different 
particles are not too close so that the inner discs centered at them do not overlap, 
a particle detection occurs when the measurement ml — m2 is high. In addition 
to cancel the low-frequency components of the single measure mi, the measure m 
yields detection thresholds independent of /i, which is valuable when the precise 
value of fji is not known. In a sense, m 2 can be considered a local estimate of /i. 
The NFA that takes into account the local context is as follows. 



Fig. 3: Measurements taken at each test location: average intensity mi in the 
measurement area (gray disc centered at (x,y)), and average intensity m 2 in the 
associated local context area (a ring around the measurement area). 


Definition 2 (NFA for the model with contrast to the context) Let u he 

an image, let T he its grid and let si and S 2 he two eireular, eoneentrie measure¬ 
ment kernels. A number of false alarms assoeiated to the measurements mi{x,y) = 
{u * Si)(x, y){i = 1,2, (x, y) ) for the naive model Ro = A/’(/i, a) is given hy 

= ( 2 ) 

where |T| is the number of points of the image grid, Af(ia;a) is the normalized 
two-dimensional Gaussian white noise, (fi = si*A/’(yn; a) and (^2 = s2^ff{ia] a) are 
the random variable following the naive model distribution in si and S 2 respeetively 
that differ in the £2 —norm, and ^c(t) = dt is the tail of the normal 

distribution. 


We say that the i — th point of the image grid T is e—meaningful, with e > 0, 
if NFA{i) < e. In the naive model, the random variable M = cfi — (1)2 follows the 
law A/’(0, (1 + (c 2 ^_i) ):^^)- In practice, for a given e that specifies the NFA and 
it is usually taken to be 1, we perform at every pixel location i the following test: 


mi — m2 
a 


> T(e, R, a). 


( 3 ) 
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where a is an estimation of the standard deviation of naive model and T(e, a) = 

er/c“^(^) 

In Equation!^ the standard deviation a is set once for all and generally has 
to be estimated. In presence of both high and low contrasted spots, this could 
be problematic since if the estimation of a is bigger than the real noise standard 
deviation, low contrasted particles would be missed. For the opposite case, since 
high intensity values have a larger dynamic range than low intensity values, as 
soon as the measure m becomes positive, it is very likely that m becomes big com¬ 
pared to the global estimate of cr. In these conditions, even in locations very bad 
centered near a real particle, the NFA is typically small leading to very widespread 
detections, while the only meaningful detection is the center of the particle. To 
solve this problem and to allow the detection, without any bias, of low and high 
contrasted particles, we relaxed the weight due to the contrast by replacing the 
global estimation of a in Eq. by a local estimation of cr in a local neighborhood 
on the candidate particle, say cr/. This choice leads to a redefinition of the naive 
models that becomes a Gaussian white noise with a local standard deviation ai. 



Radius 


(a) (b) 

Fig. 4: Left: measures of meaningfulness on a perfect Gaussian particle of cr = 
10 pixels. Right: corresponding maximum of the gradient (in black) and of the 
function (R) (in white, Ropt = 1.45cr) 


4.1 Particle spreading estimate and hiding process 

The proposed detection approach leads to a natural particle spreading estimation 
which allows to quantify the amount of fiuorescence of a particle. For a given 
particle location, once a is set, the ratio is a function of R, the inner 

radius of the model. This function depends solely on the local contrast and its 
maximum gives an estimate of the particle spreading. In fact, as it can be observed 
in Fig. 4 for a 2D-circular Gaussian particle, the maximum of the ratio is 












10 


Mariella Dimiccoli et al. 


proportional to the spread of the Gaussian. Let 


Ropt = argmax{ 
R 


mi — m2 


(R)) 


( 4 ) 


be the estimation of the particle spreading. 



Fig. 5: (a) Original image, (b) 0.1-meaningful detection, (c) 2-pass 0.1-meaningful 
detections with hiding process {R = 2.5 pixels, a = 2) 



The function {R) is drawn in Figj^ (a) for a local model centered on 



a 2D-circular Gaussian function with standard deviation a = 10. In this case 
Ropt = 1.45. On real images, we compute the particle center locations for a given 
model and a given e as the gravity center of e-meaningful connected components. 
The function (R) is computed on each detected center and its maximum 

Ropt gives a radius estimation. 

The particle spreading estimation is also useful to improve the algorithm per¬ 
formances: the detection of low contrasted particles may fail when they are very 
close to large and high contrasted particles. To overcome such a problem, we con¬ 
sider that almost all fluorescence values of a given particle (x, y) are contained in 
the 95% confldence interval, i.e inside a circle with a radius i? 2 <; = 2/1 AbRopt- 
Then the algorithm removes any detected particle inside the 95% confldence inter¬ 
val and redoes the detection (see Fig.[^. Practically, pixel values inside the circle 
C(x, y, when computing mi, m 2 and ai on close candidate particle locations, 
are set to the value of the background. Since this process deforms the pattern of 
the local model near former detections, the minimal distance between a hidden 
particle and the inner part of a new local model should be at least one pixel. 

4.2 Sub-pixel reflnement 

In fluorescence imaging the target of interest is usually smaller than the pixel size 
and all that is observed is the instrument’s PSF centered on the particle’s position 
that typically spreads several pixels wide. Therefore, in the study of particle tem¬ 
poral dynamics, it is of crucial interest to provide efficient solutions for sub-pixel 
detection. For instance, in the study of bacteria aging, the molecular components 
of interest are protein aggregates accumulated near bacteria boundaries. This par¬ 
ticular location makes very ambiguous the correspondence between aggregates and 
cells, and sub-pixel accuracy becomes crucial to disambiguate this association. To 
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noise standard deviation 


Fig. 6: Mean localization over 300 synthetic images as a function of the standard 
deviation of the Gaussian white noise. 


achieve subpixel accuracy, we propose to refine the location of each detected par¬ 
ticle, say X, by computing a weighted average of pixel coordinates in a circular 
neighborhood. The sub-pixel location, x is given by the following equation 


X = 


1 


u{y)w{y) 


( 5 ) 


where Mx is the circular neighborhood of the particle x having cardinality 
\Mx\^ w(y) = u(y) — nix if u(y) — nix > 0, w{y) = 0 otherwise, being nix is the 
median intensity value of all pixels belonging to Mx- The median value nix is an 
estimation of the local background intensity: only pixels having intensities bigger 
than the local background intensity are considered in the weighted average. This 
avoids to include into the average pixel values corresponding to very close particles. 
To validate the proposed sub-pixel refinement, we considered a set of 300 synthetic 
images with a single particle at a sub-pixel location x, simulated by a Gaussian 
bi-dimensional signal of standard deviation as and noised with a white noise of 
standard deviation an- We computed the mean error over the set of images as a 
function of the parameter as and of the radius r of the circular neighborhood Mx - 
Using a radius r = 6, we achieved a precision of 1/10 of pixel on synthetic images 
even on very poor contrasted particles (see Fig. [^. 
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5 A-contrario particle temporal linking 

the a-contrario model for particle temporal linking, first introduced in [46] as AS- 
TRE, estimates the probability of observing a trajectory in random data through a 
probabilistic criterion which combines data statistics, such as the number of images 
and particles, with trajectory characteristics such as trajectory length, particle 
density, and smoothness. This criterion is used to drive a dynamic programming 
algorithm 0 which sequentially extract the most meaningful trajectories globally 
in time, while guarantying that no trajectory will be found in random data. In 
the following we give a brief introduction to this method and to its more recent 
variant, called CUTASTRE [1]. 


5.1 Expliciting the a-contrario model for trajectory detection 

As for any othera-contrario-based method, the a-contrario approach for the detec¬ 
tion of trajectories, is grounded on two elements: the naive model and a statistical 
measurement function able to characterize the visual saliency of the trajectories. 

In the following, we assume that we are given K images /i,..., each image Ik 
containing N points Xf, corresponding to particles that have already been 

detected in each image of the sequence. We assume that the detections are noised 
so that the observed points may correspond to spurious particles and that some 
particles could have been missed. Intuitively, by Helmholtz principle, we should 
not see trajectories appearing in the realizations of the naive model. In the case 
of trajectory detection, the naive model is an uniform and iid draw of N points 
in each of the K frames. Before introducing the statistical measurement function, 
we define the structures of interest. 

Definition 3 (Trajectories without holes) A trajeetory T of length i start¬ 
ing at frame ko is a tuple T = (ko,ii^ --Ai), where 1 < ip < N for all p and 
1 < i < K — ko 1. We will denote by T the set of all trajeetories. There 
is a natural equivalenee between a trajeetory T G T and the tuple of variables 
Xt = {X^f ^ that we shall therefore sometimes abusively eall a (ran¬ 

dom) trajeetory too. 

The statistical measurement function associated to a trajectory T = (X^f , ^ 

with length /, where Xf is the i-th point of frame /c, is its maximal acceleration 
amax(t) = max 3 <^<£_i ||yz+i — 2yi +yi_i||. The amount of surprise when observ¬ 
ing a trajectory T of length / and acceleration 6 := a{t) is estimated by the upper 
bound 

P«„(a(T0<<5)<(7r5V|r2|)'-2 (6) 

where Q is the image domain. The NEA can than be computed thanks to Lemma 
1 in [29] . 

Trajectories are extracted iteratively. At each iteration, the value of the tra¬ 
jectory with minimal NEA among all trajectories, say m, is computed through a 
dynamic programming strategy and compared to e, the maximal NEA value of a 
trajectory that the user wants to extract (usually one chooses e = 1). If m < e, all 
points corresponding to the trajectory are removed from the sequence. This process 
is repeated until no trajectory with NEA less than e can be found anymore. The 
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ASTRE algorithm is global in time and has an unique parameter e. However, it 
has quadratic time and memory complexity with respect to the number of frames, 
that may be prohibitive for some applications involving long image sequences (say 
more than 1000 frames). Abergel and Moisan [1] proposed a modified version of 
ASTRE, called CUTASTRE with 0{K) complexity, which preserves the rigorous 
control of false detections in pure noise offered by ASTRE achieving very similar 
performances. This important complexity reduction is obtained at the cost of two 
additional parameters: the temporal chunk and overlap sizes, that in general, as 
proved in [T], are easy to be set since their are related to the smoothness of the 
trajectory. In addition, CUTASTRE is not currently able to deal with trajectories 
with holes as ASTRE does. 


6 Experimental Results 

6.1 Comparative evaluations of particle detection 

To evaluate quantitatively the performances of the proposeda-contrario particle 
detection method, we used the same experimental setup proposed by Smal et 
al. [H]. In this work, seven unsupervised (WMP, MSVST, TH, MTH, SEE, HD, 
IDE) and two supervised methods (AB and EDA) were compared in terms of true 
positive and false positive rate, taking into account also the methods’s sensitivity 
to parameter changes and data quality. 


Table 1: Optimal NEA parameter and corresponding performances of the a- 
contrario particle detection for different SNR 


Image Type 

NFA 

TPR 

FPR 

SNR = 4 

A 

0.3 

1 

0.005 

B 

0.3 

0.96 

0.036 

C 

0.1 

1 

0.005 

SNR = 3 

A 

0.3 

1 

0.009 

B 

0.3 

0.98 

0.031 

C 

0.1 

1 

0.006 

SNR = 2 

A 

0.3 

1 

0.008 

B 

0.3 

0.99 

0.028 

C 

0.1 

0.97 

0.005 


Table 2: Radii parameters used for each type of image for all SNR 


Image Type 

R 

OL 

A 

3 

2 

B 

2 

3.5 

C 

3 

1.25 
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6AA Simulated image dataset and performanee measures 

We used the publicly available Imaged plugin syndata.jar to generate three types 
of image sets, each for SNR ranging from 2 to 4. Each type of image set corre¬ 
sponds to an uniform background (type A), a gradient background (type B) and a 
non-uniform background (type C). The performances of the a-contrario detection 
methods were evaluated by accumulating the numbers of true positive (TP) and 
false negative (FN) for 16 images, each containing 256 ground truth objects and 
averaging the results over all objects. As in m, the tolerance was fixed to 4 pix¬ 
els. To compare the algorithms were used two main measures: 1) the true-positive 
rate (TPR) defined as TPR = NTP/ {NTP + NFN), where NTP stands for 
number of true positives and NFN is the number of false negatives defined as 
NFN = NO — NTP, being A^O the number of objects in the ground truth; 2) 
the modified false positive rate defined as FPR* = NFP/{NTP + NFN), where 
NFP IS, the number of false positives. The value of TPR for the optimal param¬ 
eters is denoted as TPR*. In addition to detection performances with optimal 
parameters, the sensitivity to parameter changes and data quality is also consid¬ 
ered. The sensitivity to parameter changes is evaluated by plotting the FROG 
curves obtained at a fixed SNR for two different values of a given parameter and 
for each image type. The sensitivity to data quality is quantified by considering 
the variation of TPR and FPR for different values of the SNR, keeping fixed the 
algorithm’s parameters for each image type. 

6.1.2 Diseussion 

In Table are shown the performances of the proposed method for the optimal 
NFA parameter for varying SNRs. These results have been obtained using fixed 
values of the internal radius R and a for each image type, which are reported in 
Table These experiments demonstrate that the proposed method is robust to 
variations of data quality. 

To evaluate the sensitivity of the proposed method to parameter changes, in 
Figj^(a) we report the FROG curves obtained at SNR = 2 for two different values 
of the internal radius R for each image type, varying a. As it can be observed, 
the performances are quite stable with respect to variations of these parameters 
and performances do not go below the 96% TPR. In Table we compare our 
results to those obtained by using FDA and ADABOOST that, in the comparison 
work of Smal et al. EH, reported the highest TPR* and the lowest sensitivity to 
parameter changes and data quality on our same dataset with respect to seven 
unsupervised methods (TH, SEF, WMP, MSVTV, MTH, HD, IDF) . The sensi¬ 
tivity of the measures TPR and FPR* to a parameter, say Id (that correspond 
to the NFA for our method and to the threshold on the size of the clusters for 
FDA and ADABOOST) is measured through the values St = —dTPRldld and 
Sf = —{dFPR*/did) at the value of the parameter for which the FPR* = 0.01 
(only 1% false positives) hereafter called optimal parameter. As it can be observed, 
the sensitivity of our method to variations of the NFA is very small and order of 
magnitude smaller than the one of AB and FDA. 

The proposed method has been extensively tested on biological images to detect 
protein aggregates in growing bacteria cultures as reported in m- In this work, 
to deal with the presence of particles having variable size, we used a multiscale 
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approach that in turn consists in using two different sizes for the inner radius. In 
such experiments, the NFA was fixed to 1, a = 2 and inner radii R varying from 2 
to 3 pixels. The same process was performed on four different growing sequences 
of about 90 images comprising 1 or 2 particles at the beginning and up to about 
150 particles at the end. In mean 15 particles per sequence were recovered thanks 
to the hiding process. 
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Fig. 7: (a) FROG curves for the a-contrario method depending on the value of the 
R and aR at SNR = 2 and optimal value of the NFA. (b) FROG curves for AB 
method depending on the value of the Vd threshold at SNR = 2 and with number 
of Haar-like features Nab = 50. (c) FROG curves for FDA method depending on 
the value of the Vd threshold at SNR = 2. Image source for (b) and (c) [61] . 
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ACONTRARIO 

Image type 

NFA 

TPR 

St 

Sp 

A 

0.3 

1 


10“^ 

B 

0.3 

0.99 

lO-'^ 

(M 

1 

o 

1—1 

C 

0.1 

0.97 

10-5 

(M 

1 

o 

1—1 

ADABOOST 

Image type 


TPR 

St 

Sf 

A 

3 

0.99 

10-^ 

10-3 

B 

31 

0.94 

.01 

10-3 

C 

30 

0.94 

.01 

10-3 

FDA 

Image type 


TPR 

St 

Sp 

A 

4.6 

0.99 

10-5 

.01 

B 

8.8 

0.99 

10-3 

.01 

C 

9.8 

0.96 

(M 

1 

o 

1—1 

.01 


Table 3: Optimal parameters for each image type at SNR = 2. is the threshold 
on the size of the clusters used by ADABOOST and FDA; NFA is the threshold 
used by the proposed method. 


From the above experiments, we can conclude that the proposed approach 
overall outperforms both supervised and unsupervised methods considered in m 
in terms of true positive ratio with only 1% false positives at very low SNR and 
in terms of robustness with respect to data quality and parameter changes. In 
addition, the proposed method has the advantage of not requiring a cumbersome 
training stage on similar data, which requires additional parameters to be esti¬ 
mated such as the number of Haar-like features. Finally, it has the advantage of 
allowing the user to directly set the parameter e, which represent a bound of the 
average number of detections that would be made in pure noise data (that is on 
spurious detections). Due to its robustness to parameter changes and to poor data 
quality, as well as to its unsupervised nature, our method is very suited to be used 
by biologists. In the next section we further validate the suitability of the proposed 
method as input for a particle temporal linking algorithm. 


6.2 Comparative evaluations of particle detection followed by temporal linking 

In this section, we evaluate the performances of the a-contrario temporal link¬ 
ing approach called CUTASTRE introduced in [1] when particles are detected 
by using the a-contrario approach proposed in this paper. Our goal is twofold: 
first, to validate the suitability of the proposed particle detection method as input 
for the task of linking particles in successive frames; second, to demonstrate the 
advantages of the a-contrario appraoch for both particle detection and temporal 
linking. For that, we used the baseline issue of the 2012 Particle Tracking Chal¬ 
lenge data (see http://www.bioimageanalysis.org/track) to which participated 14 
teams. Each team applied his method independently on a common dataset and 
evaluated the results using a set of commonly evaluation criteria. 
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6.2.1 Simulated image data sets 

Since the ground truth is generally not available for real biological data and man¬ 
ual annotation by human observers is subjective, costly and prone to bias m, 
the authors simulated image data for the challenge. The simulated data, together 
with their corresponding ground truth, take into account three factors that usu¬ 
ally have large influence on tracking results: the particle dynamics characterizing 
a biological scenario, the particle density in a flxed held of view and the particle 
signal relative to noise. Four biological scenarios were simulated, including near¬ 
circular particles showing a Brownian motion (VESICLES), near-circular particles 
switching between Brownian and directed motion models (RECEPTORS), near¬ 
circular particle captured in 3D switching between Brownian and directed motion 
models (VIRUS), elongated particles showing near-constant velocity motion (MI- 
CROTUBLES). Eor each biological scenario three levels of particle density were 
considered (low 100 particles, mid 500 particles, high 1000 particles), and four 
SNRs (1,2,4,7). Being the proposed algorithms suitable only for 2D detection and 
tracking of symmetric Gaussians, we considered only 2 of the 4 aforementioned 
scenarios, VESICLE and RECEPTOR, amounting to a data set consisting of 24 
image sequences. It is worth to mention that the 3D scenario represented by VIRUS 
is the easiest from a temporal-linking point of view since these particles shows a 
near-constant velocity motion and VIRUS show the same motion mode than RE¬ 
CEPTORS, being the only difference the shape of the particles. 

6.2.2 Performanee measures 

A set of flve complementary criteria that give a complete and intuitive characteri¬ 
zation of the tracking results when tracking with a varying number of particles in 
a cluttered environment were used to evaluate the tracking performances. In the 
following, we give an intuitive explanation of each of them and we refer the reader 
to the work of nHi for further details. Let us deflne a track 9 that exists from time 
tstart to time tend as a temporal series of subsequent spatial positions, say the set 
0 = {9{t) = {x{t), with t = tstart, ■■■, tend- 

The distance between two tracks is deflned as 

T-l 

a{9l,92)= V - feft)||2,e 

t=0 

where T is the lenght of the image sequence and the distance between two positions 
is deflned as 


|| 6 »i(i) - 6 » 2 (t)|| 2 ,e = min(||6»i(t) - 6 » 2 (t)|| 2 ,e) 

where || - || is the standard £2 norm of TZ^ and e G Ti+. This measure limits the 
penalization for tracks that are more than e apart to e. The parameter e was flxed 
to flve in the particle tracking challenge. 

The distance between two track sets is deflned as follows. Let X be the 
set of ground-truth tracks and, Y the estimated set of tracks and Y the extended 
version of Y with dummy tracks. Denoting by i? the set of tracks that can be 
obtained by taking |X| elements from V, the distance between X and an element 
Z from i? is deflned as the sum of the distances between |X| pairs given by the 
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ordering of the two sets. The distance between X and Y is then defined as the 
minimum distance between X and all possible Z: 

1^1 

d(X,Y) = min 

^ k=l 

By using the above definitions of distance, five metrics used for evaluation are 
defined as follows: 

1. The normalized score: a(^ 1 , ^ 2 ) = l — ^ 2 )/d(^i, where ^ denotes the set 

of 1^11 dummy trucks. It measures the overall degree of matching of groundtruth 
and estimated tracks without taking into account spurious tracks. 

2. The criterion 13{ 61 ^ 62 ) = S) + measures the overall degree of 

matching of groundtruth and estimated tracks with a penalization of nonpaired 
estimated tracks. 

3. Jaccardi similarity index for positions JSC = tp+fn+fp • inter¬ 

val [0,1] and takes value 1 only if there are not spurious tracks in Y and all 
positions pairs in (X, Z*) are matching. 

4. Jaccard similarity coefficient for entire tracks JSCe = TPe+^N^e+FPe ' 

in the interval [0,1] and takes value 1 only if there are not spurious tracks in 
Y and Z* does not contain dummy tracks. 

5. Root Mean Squared Error (RMSE) that indicates the overall localization accu¬ 
racy of matching points in the optimally paired tracks by using the Euclidean 
distance 

6.2.3 Discussion 

In this section we discuss comparative performances of the a-contrario approach 
against 14 methods that partecipated to 2012 Particle Tracking Challenge. Since in 
the dataset only trajectories without holes are considered, we used CUTASTRE [1] 
which cannot handle trajectories with holes but has linear complexity. Eor particle 
detection we used a multiscale approach consisting of using different values of the 
internal radii (from 1.5 to 3). We varied the algorithm parameters in a small range 
and kept the best result for each scenario. 

The methods for particle detection used in the challenge can be roughly grouped 
in four classes (see Table |^: 1) thresholding or local maxima selection (meth¬ 
ods 2,3,4,9); 2) linear (Gaussian, Laplacian of Gaussian and difference of Gaus¬ 
sian) and nonlinear model fitting including morphological processing (methods 
6,7,11,12,13,14) 4) centroid estimation scheme (method 1) or a combination of 
them (methods 5,8,10). The methods for temporal linking include deterministic 
(method 1, 6, 8, 9, 12, 13,) and probabilistic approaches (methods 2,3,4,5,7,10,11,14) 
and they were introduced in section In Eigj^ Fig(^ and EigfTol comparisons 
are shown for low, mid and high particle density scenarios respectively. The left 
and the right columns show the performances in terms of the five measures defined 
above for REGEPTORS and VESGIGLES respectively as a function of the SNR. 
As it can be observed, the performaces of all methods, within a given scenario, 
depends on particle density and SNR. Analyzing trends, it emerges that generally 
the performances of the a-contrario approach, called Method 15 in the legend, de¬ 
creases less strongly when the SNR goes from 4 to 2 with respect other methods. 






20 


Mariella Dimiccoli et al. 


The metrics a and (3 differ slighly independently on the particle density and on 
the kind of particle motion (Brownian for VESICLE and Brownian switching to 
directed motion models for RECEPTORS), meaning that the a-contrario method 
gives a low number of spuriuos tracks. As general trend, these measures for the 
a-contrario approach decrease while increasing particle density, even if they are 
above the state of the art. This can be understood considering that the parti¬ 
cle detection algorithm assumes that the minimal distance between two particles 
is at least of one pixel, an assumption that often does not hold in high density 
scenario. The Jaccardi similarity index for positions is clearly showing the best 
performances even for high particle density and low SNRs, for both VESICLE 
and RECEPTORS. This can mainly be due to the low number of spurious tracks, 
which are penalized by this measure, indicating that the tracking algorithm is ro¬ 
bust to The Jaccard similarity coefficient for entire tracks has in general higher 
value for RECEPTORS than for VESICLES with respect to other methods. This 
means that the a-contrario particle temporal linking is able to cope with the more 
complex dynamics of these particles better than competing methods. The local¬ 
ization accuracy, expressed in terms of RMSE is comparable with Method 12 of 
the state of the art which uses parabolic fitting for particle detection. However, 
it should be taken into account that this measure is an average over all corrected 
matching pairs, which are not too much for this method, as demonstrated by its 
performances in terms of a. In general, the better performances of the a-contrario 
approach are mainly due to its control on the number of false alarms, which are 
penalized in 3 out of 5 the evaluation criteria. 


Table 4: Detection and temporal linking methods compared in m 
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Met. 

Detection approach 

Linking approach 

Refs. 

1 

Iterative intensity- 

weighted centroid 

calculation 

Gombinatorial tracking 
(deterministic) 

m 

2 

Adapt at ive local- 

maxima selection 

Multiple hypothesis 

tracking (probabilistic) 

El HO] 

3 

Maxima after thresh¬ 
olding two-scale 

wavelet products 

Multiple hypothesis 

tracking (Probabilistic) 

HZllllllS] 

4 

Adaptive Otsu Thresh¬ 
olding 

Multitemporal associa¬ 
tion tracking 

EiEH] 

5 

Thresholding + cen¬ 
troid calculation 

Kalman filtering + prob¬ 
abilistic data association 

EIES] 

6 

Lorentzian function fit¬ 
ting to structures above 
noise level 

Dynamic programming 

m 

7 

Gaussian mixture 

model fitting 

Multiple hypothesis 

tracking 

m 

8 

Watershed-based 
clump splitting and 
parabola fitting 

Viterbi algorithm on 
state-space representa¬ 
tion 

Haul] 

9 

Maxima with pixel pre¬ 
cision 

Nearest neighbor + 
global optimization 

Eina 

10 

Histogram-based 
thresholding and 

Gaussian fitting 

Gaussian template 

matching 

[Hnipnip] 

11 

Gaussian fitting 

Sequential multiframe 
assignment 

gSlIMlISS] 

12 

Parabolic fitting to lo¬ 
calized maxima 

Linear assignment prob¬ 
lem 

eies 

13 

Watershed-based 
clump splitting 

Nearest neighbor 

Eina 

14 

Morphological opening- 
based clump splitting 

Nearest neighbor + 
Kalman filtering 

EiEH] 
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Fig. 8: Performance measures for low-density particle scenarios: (a) RECEPTORS, 
(b) VESICLES. Erom up to down are show the metrics a, (3^ JSC, JSCt and 
RMSE. Method 15 is the proposed approach. 


Method 1 
— Method 2 
Method 3 
Method 4 
Method 5 
Method 6 
Method 7 
Method 8 
—^ Method 9 
Method 10 
Method 11 
— Method 12 
Method 13 
Method 14 
Method 15 


- Method 1 

- Method 2 

- Method 3 

- Method 4 

■ Method 5 

- Method 6 

■ Method 7 

- Method 8 

- Method 9 

- Method 10 

- Method 11 

- Method 12 

- Method 13 

- Method 14 

- Method 15 

































































































































Title Suppressed Due to Excessive Length 


23 
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Fig. 9: Performance measures for mid-density particle scenarios: (a) RECEPTORS, 
(b) VESICLES. Erom up to down are show the metrics a, /3, JSC^ JSCt and 
RMSE. Method 15 is the proposed approach. 
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Fig. 10: Performance measures for high-density particle scenarios: (a) RECEP¬ 
TORS, (b) VESICLES. Erom up to down are show the metrics a, (3^ JSC, JSCt 
and RMSE. Method 15 is the proposed approach. 
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7 Conclusions 

This paper has shown the advantages of the a-contrario framework for the spatial 
detection and tracking of near-circular particles in fluorescent time-lapse images. 
Comparative evaluations under different biological scenarios and varying experi¬ 
mental conditions, both of the particle detection method and of particle detection 
followed by temporal linking, have shown that the proposed approach outperforms 
the state-of-the-art in terms well established performance measures. In addition 
to better performances for very low SNR, the a-contrario approach provides three 
additional advantages: a rigorous control of false detections in pure noise, which is 
important to avoid the corruption of quantitative analysis in biological data; low 
sensitivity to parameters changes; no need of a costly training stage that, could 
possibly introduce a bias. These characterisitcs make the proposed algorithms par¬ 
ticularly suited to be used by biologists. The Imaged plugin of the particle detection 
algorithm can be found online at http://fluobactracker.inrialpes.fr/, Future work 
will extend this framework to handle 3D data and elongated particles. 
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